#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TEllipse.h"
#include "TClonesArray.h"

using namespace std;

void readtest(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  int ene =100;

  TFile *fread=new TFile("../outputFiles/test.root");
  //  TFile *fread=new TFile(Form("../outputFiles/Ar_0V_%dkeV.root",ene));
  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  int nPhoton;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("nPhoton",&nPhoton);

  double LowEdge = -500;
  double HighEdge = 100000;

  TH1F *hJPMTx=new TH1F("hJPMTx","hJPMTx",100,-150,150);
  TH1F *hJPMTy=new TH1F("hJPMTy","hJPMTy",100,-150,150);
  TH1F *hJPMTz=new TH1F("hJPMTz","hJPMTz",100,-650,-550);
  TH1F *hMPMTx=new TH1F("hMPMTx","hMPMTx",100,-150,150);
  TH1F *hMPMTy=new TH1F("hMPMTy","hMPMTy",100,-150,150);
  TH1F *hMPMTz=new TH1F("hMPMTz","hMPMTz",100,-650,-550);

  const int NEV=400;
  TH1F *hJPMTt[NEV];
  TH1F *hMPMTt[NEV];
  TH2F *hPMTxy[NEV];

  TH2F *hPosHit_All=new TH2F("hPosEff_All","hPosEff_All",30,-150,150,35,-650,-300);
  TH2F *hPosHit_J=new TH2F("hPosEff_J","hPosEff_J",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M1=new TH2F("hPosEff_M1","hPosEff_M1",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M2=new TH2F("hPosEff_M2","hPosEff_M2",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M3=new TH2F("hPosEff_M3","hPosEff_M3",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M4=new TH2F("hPosEff_M4","hPosEff_M4",30,-150,150,35,-650,-300);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;
  if(nev>NEV) nev=NEV;

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    tree->GetEntry(iev);

    hJPMTt[iev]=new TH1F(Form("hJPMTt_%d",iev),Form("hJPMTt_%d",iev),(HighEdge-LowEdge)/10,LowEdge,HighEdge);
    hMPMTt[iev]=new TH1F(Form("hMPMTt_%d",iev),Form("hMPMTt_%d",iev),(HighEdge-LowEdge)/10,LowEdge,HighEdge);
    hPMTxy[iev]=new TH2F(Form("hPMTxy_%d",iev),Form("hPMTxy_%d",iev),100,-150,150,100,-150,150);

    int njhit = JPMTa->GetEntriesFast();
    for(int ijhit=0;ijhit<njhit;ijhit++){
      PMT *jpmt = (PMT*)JPMTa->At(ijhit);
      hJPMTx->Fill( jpmt->xpos );
      hJPMTy->Fill( jpmt->ypos );
      hJPMTz->Fill( jpmt->zpos );
      hJPMTt[iev]->Fill( jpmt->time );
      hPMTxy[iev]->Fill( jpmt->xpos,jpmt->ypos);
    }
    int nmhit1 = MPMT1a->GetEntriesFast();
    for(int imhit=0;imhit<nmhit1;imhit++){
      PMT *mpmt1 = (PMT*)MPMT1a->At(imhit);
      hMPMTx->Fill( mpmt1->xpos );
      hMPMTy->Fill( mpmt1->ypos );
      hMPMTz->Fill( mpmt1->zpos );
      hMPMTt[iev]->Fill( mpmt1->time );
      hPMTxy[iev]->Fill( mpmt1->xpos,mpmt1->ypos);
    }

    int nmhit2 = MPMT2a->GetEntriesFast();
    int nmhit3 = MPMT3a->GetEntriesFast();
    int nmhit4 = MPMT4a->GetEntriesFast();


    int iz = iev / 21 % 18;
    int ir = iev % 21;
    float zpos = iz*10.-545;
    float rpos = ir*10.-120;

    hPosHit_All->Fill( rpos,zpos,(njhit+nmhit1+nmhit2+nmhit3+nmhit4) );
    hPosHit_J->Fill( rpos,zpos,njhit );
    hPosHit_M1->Fill( rpos,zpos,nmhit1 );
    hPosHit_M2->Fill( rpos,zpos,nmhit2 );
    hPosHit_M3->Fill( rpos,zpos,nmhit3 );
    hPosHit_M4->Fill( rpos,zpos,nmhit4 );

  }

  /*
  TCanvas *c1=new TCanvas("c1","c1",1000,800);
  c1->Divide(3,2);
  c1->cd(1);
  hJPMTx->Draw();
  c1->cd(2);
  hJPMTy->Draw();
  c1->cd(3);
  hJPMTz->Draw();
  c1->cd(4);
  hMPMTx->Draw();
  c1->cd(5);
  hMPMTy->Draw();
  c1->cd(6);
  hMPMTz->Draw();

  TCanvas *c2=new TCanvas("c2","c2",1200,1000);
  c2->Divide(5,4);
  TEllipse *eJPMT=new TEllipse(0.,0.,32.,32.);
  TEllipse *eMPMT=new TEllipse(-2*40,0.,23.,23.);
  TCanvas *c3=new TCanvas("c3","c3",1200,1000);
  c3->Divide(5,4);

  for(int i=0;i<20;i++){
    c2->cd(i+1);
    hPMTxy[i]->Draw("colz");
    eJPMT->SetFillStyle(0);
    eMPMT->SetFillStyle(0);
    eJPMT->Draw();
    eMPMT->Draw();
    c3->cd(i+1);
    hJPMTt[i]->SetXTitle("time (ns)");
    hJPMTt[i]->SetYTitle("#P.E.");
    hJPMTt[i]->Draw();
  }
  */

  gStyle->SetOptStat(000000);

  TCanvas *c4=new TCanvas("c4","C4",1200,900);
  c4->Divide(3,3);
  c4->cd(5);
  hPosHit_J->Draw("colz");
  c4->cd(4);
  hPosHit_M1->Draw("colz");
  c4->cd(8);
  hPosHit_M2->Draw("colz");
  c4->cd(6);
  hPosHit_M3->Draw("colz");
  c4->cd(2);
  hPosHit_M4->Draw("colz");



}

void readtest2(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

    TFile *fread=new TFile("../outputFiles/test.root");
  //  TFile *fread=new TFile("../outputFiles/10000photon.root");
  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  int nPhoton;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("nPhoton",&nPhoton);

  double LowEdge = -500;
  double HighEdge = 100000;

  TH1F *hJPMTx=new TH1F("hJPMTx","hJPMTx",100,-150,150);
  TH1F *hJPMTy=new TH1F("hJPMTy","hJPMTy",100,-150,150);
  TH1F *hJPMTz=new TH1F("hJPMTz","hJPMTz",100,-650,-550);
  TH1F *hMPMTx=new TH1F("hMPMTx","hMPMTx",100,-150,150);
  TH1F *hMPMTy=new TH1F("hMPMTy","hMPMTy",100,-150,150);
  TH1F *hMPMTz=new TH1F("hMPMTz","hMPMTz",100,-650,-550);

  const int NEV=400;
  TH1F *hJPMTt[NEV];
  TH1F *hMPMTt[NEV];
  TH2F *hPMTxy[NEV];

  TH2F *hPosHit_All=new TH2F("hPosEff_All","hPosEff_All",30,-150,150,35,-650,-300);
  TH2F *hPosHit_J=new TH2F("hPosEff_J","hPosEff_J",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M1=new TH2F("hPosEff_M1","hPosEff_M1",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M2=new TH2F("hPosEff_M2","hPosEff_M2",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M3=new TH2F("hPosEff_M3","hPosEff_M3",30,-150,150,35,-650,-300);
  TH2F *hPosHit_M4=new TH2F("hPosEff_M4","hPosEff_M4",30,-150,150,35,-650,-300);
  TH1F *hMeanAll=new TH1F("hMeanAll","hMeanAll",1000,0,1);
  TH1F *hMeanJ=new TH1F("hMeanJ","hMeanJ",1000,0,1);
  TH1F *hMeanM1=new TH1F("hMeanM1","hMeanM1",1000,0,1);
  TH1F *hMeanM2=new TH1F("hMeanM2","hMeanM2",1000,0,1);
  TH1F *hMeanM3=new TH1F("hMeanM3","hMeanM3",1000,0,1);
  TH1F *hMeanM4=new TH1F("hMeanM4","hMeanM4",1000,0,1);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;
  //if(nev>NEV) nev=NEV;

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    tree->GetEntry(iev);

    hJPMTt[iev]=new TH1F(Form("hJPMTt_%d",iev),Form("hJPMTt_%d",iev),(HighEdge-LowEdge)/10,LowEdge,HighEdge);
    hMPMTt[iev]=new TH1F(Form("hMPMTt_%d",iev),Form("hMPMTt_%d",iev),(HighEdge-LowEdge)/10,LowEdge,HighEdge);
    hPMTxy[iev]=new TH2F(Form("hPMTxy_%d",iev),Form("hPMTxy_%d",iev),100,-150,150,100,-150,150);

    int njhit = JPMTa->GetEntriesFast();
    for(int ijhit=0;ijhit<njhit;ijhit++){
      PMT *jpmt = (PMT*)JPMTa->At(ijhit);
      hJPMTx->Fill( jpmt->xpos );
      hJPMTy->Fill( jpmt->ypos );
      hJPMTz->Fill( jpmt->zpos );
      hJPMTt[iev]->Fill( jpmt->time );
      hPMTxy[iev]->Fill( jpmt->xpos,jpmt->ypos);
    }
    int nmhit1 = MPMT1a->GetEntriesFast();
    for(int imhit=0;imhit<nmhit1;imhit++){
      PMT *mpmt1 = (PMT*)MPMT1a->At(imhit);
      hMPMTx->Fill( mpmt1->xpos );
      hMPMTy->Fill( mpmt1->ypos );
      hMPMTz->Fill( mpmt1->zpos );
      hMPMTt[iev]->Fill( mpmt1->time );
      hPMTxy[iev]->Fill( mpmt1->xpos,mpmt1->ypos);
    }
    int nmhit2 = MPMT2a->GetEntriesFast();
    int nmhit3 = MPMT3a->GetEntriesFast();
    int nmhit4 = MPMT4a->GetEntriesFast();

    int iz = iev / 21 % 18;
    int ir = iev % 21;
    float zpos = iz*10.-545;
    float rpos = ir*10.-120;

    hPosHit_All->Fill( rpos,zpos,(njhit+nmhit1+nmhit2+nmhit3+nmhit4)/(float)nPhoton );
    hPosHit_J->Fill( rpos,zpos,njhit/(float)nPhoton );
    hPosHit_M1->Fill( rpos,zpos,nmhit1/(float)nPhoton );
    hPosHit_M2->Fill( rpos,zpos,nmhit2/(float)nPhoton );
    hPosHit_M3->Fill( rpos,zpos,nmhit3/(float)nPhoton );
    hPosHit_M4->Fill( rpos,zpos,nmhit4/(float)nPhoton );

    hMeanAll->Fill( (njhit+nmhit1+nmhit2+nmhit3+nmhit4)/(float)nPhoton );
    hMeanJ->Fill( njhit/(float)nPhoton );
    hMeanM1->Fill( nmhit1/(float)nPhoton );
    hMeanM2->Fill( nmhit2/(float)nPhoton );
    hMeanM3->Fill( nmhit3/(float)nPhoton );
    hMeanM4->Fill( nmhit4/(float)nPhoton );
    
  }

  cout<<" All Eff(Mean) = "<<100*hMeanAll->GetMean()<<" %"<<endl;
  cout<<"JPMT Eff(Mean) = "<<100*hMeanJ->GetMean()<<" %"<<endl;
  cout<<"MEG1 Eff(Mean) = "<<100*hMeanM1->GetMean()<<" %"<<endl;
  cout<<"MEG2 Eff(Mean) = "<<100*hMeanM2->GetMean()<<" %"<<endl;
  cout<<"MEG3 Eff(Mean) = "<<100*hMeanM3->GetMean()<<" %"<<endl;
  cout<<"MEG4 Eff(Mean) = "<<100*hMeanM4->GetMean()<<" %"<<endl;

/*
  TCanvas *c1=new TCanvas("c1","c1",1000,800);
  c1->Divide(3,2);
  c1->cd(1);
  hJPMTx->Draw();
  c1->cd(2);
  hJPMTy->Draw();
  c1->cd(3);
  hJPMTz->Draw();
  c1->cd(4);
  hMPMTx->Draw();
  c1->cd(5);
  hMPMTy->Draw();
  c1->cd(6);
  hMPMTz->Draw();

  TCanvas *c2=new TCanvas("c2","c2",1200,1000);
  c2->Divide(5,4);
  TEllipse *eJPMT=new TEllipse(0.,0.,32.,32.);
  TEllipse *eMPMT=new TEllipse(-2*40,0.,23.,23.);
  TCanvas *c3=new TCanvas("c3","c3",1200,1000);
  c3->Divide(5,4);

  for(int i=0;i<20;i++){
    c2->cd(i+1);
    hPMTxy[i]->Draw("colz");
    eJPMT->SetFillStyle(0);
    eMPMT->SetFillStyle(0);
    eJPMT->Draw();
    eMPMT->Draw();
    c3->cd(i+1);
    hJPMTt[i]->SetXTitle("time (ns)");
    hJPMTt[i]->SetYTitle("#P.E.");
    hJPMTt[i]->Draw();
  }
  */
  gStyle->SetOptStat(000000);

  TCanvas *c4=new TCanvas("c4","C4",1200,900);
  c4->Divide(3,3);
  c4->cd(5);
  hPosHit_J->Draw("colz");
  c4->cd(4);
  hPosHit_M1->Draw("colz");
  c4->cd(8);
  hPosHit_M2->Draw("colz");
  c4->cd(6);
  hPosHit_M3->Draw("colz");
  c4->cd(2);
  hPosHit_M4->Draw("colz");

}


void readtest3(){

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  TFile *fread=new TFile("../outputFiles/test.root");
  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  TClonesArray *Particlea=new TClonesArray("Particle");
  int nPhoton;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("Particleb",&Particlea);
  tree->SetBranchAddress("nPhoton",&nPhoton);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  TGraph *XYplane=new TGraph();
  TGraph *XZplane=new TGraph();

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    Particlea->Clear();
    tree->GetEntry(iev);

    int nstep = Particlea->GetEntriesFast();
    for(int istep=0;istep<nstep;istep++){
      Particle *particle1 = (Particle*)Particlea->At(istep);
      if(iev==2){
	if(istep==0){
	  XYplane->SetPoint(istep,particle1->xpre,particle1->ypre);
	  XZplane->SetPoint(istep,particle1->xpre,particle1->zpre);
	}
	XYplane->SetPoint(istep+1,particle1->xpos,particle1->ypos);
	XZplane->SetPoint(istep+1,particle1->xpos,particle1->zpos);
	cout<<istep<<" "<<particle1->particlename<<" "<<particle1->xpos<<" "<<particle1->ypos<<" "<<particle1->zpos<<endl;
      }
    }

  }

  TEllipse *eInner=new TEllipse(0.,0.,150.,150.);
  eInner->SetFillStyle(0);
  TLine *lFid1=new TLine(-133,-360,-133,-550);
  TLine *lFid2=new TLine(-133,-550,  94,-550);
  TLine *lFid3=new TLine(  94,-550,  94,-360);
  TLine *lFid4=new TLine(  94,-360,-133,-360);

  
  TCanvas *c1=new TCanvas("c1","c1",1200,600);
  c1->SetGrid();
  c1->Divide(2,1);
  c1->cd(1);
  TH1 *frame1=gPad->DrawFrame(-150,-150,150,150);
  frame1->Draw();
  XYplane->SetLineColor(4);
  XYplane->SetMarkerColor(4);
  XYplane->SetMarkerStyle(21);
  XYplane->SetMarkerSize(1);
  XYplane->Draw("pl");
  eInner->Draw();
  
  c1->cd(2);
  TH1 *frame2=gPad->DrawFrame(-150,-650,150,-300);
  frame2->Draw();
  XZplane->SetLineColor(4);
  XZplane->SetMarkerColor(4);
  XZplane->SetMarkerStyle(21);
  XZplane->SetMarkerSize(1);
  XZplane->Draw("pl");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();


}
